Reservoir simulation is a process of inferring the behavior of a real reservoir from the performance of a model of that reservoir. Because mass transfer and fluid flow processes in petroleum reservoirs are so complex, reservoir simulations can only be done using computers. Computer programs that perform calculations to simulate reservoirs are called reservoir simulators. The objective of reservoir simulation is to understand the complex chemical, physical, and fluid flow processes occurring in a petroleum reservoir sufficiently well to be able to predict future behavior of a reservoir and to maximize recovery of hydrocarbons. The reservoir simulator can solve reservoir problems that are not solvable in any other way. For example, a reservoir simulator can predict the consequences of reservoir management decisions. Reservoir simulation often refers to the hydrodynamics of flow within a reservoir, but in a larger sense it also refers to the total petroleum system which includes the reservoir, the surface facilities, and any interrelated significant activity.
FIG. 1 illustrates schematically four basic steps in conventional reservoir simulation of a petroleum reservoir. The first step (step 1 in FIG. 1) is to construct a mathematical model of a real reservoir based on the chemical, physical, and fluid flow processes occurring in the reservoir. The mathematical model will comprise a set of nonlinear partial differential equations. The second step (step 2 in FIG. 1) involves discretization of the reservoir in both time and space. Space is discretized by dividing the reservoir into suitable cells with each cell having a set of nonlinear finite difference equations. The third step (step 3 in FIG. 1) is to linearize the nonlinear terms that appear in the nonlinear finite difference equations and, based on this linearization, construct linear algebraic equations. The fourth step (step 4 in FIG. 1) is to solve the linear algebraic equations. The solution provides a prediction of reservoir behavior, which enables a petroleum engineer to predict reservoir performance, including the rate at which the reservoir can be produced. The accuracy of the model can be checked against the history of the reservoir after the model has been subjected to a simulated recovery process.
In developing a mathematical model (step 1 in FIG. 1), the following principles apply:
(1) Anything that enters or leaves the reservoir must cross a boundary. PA1 (2) At some initial time the reservoir can be described by some set of conditions. PA1 (3) The processes that occur within the reservoir obey known physical laws and consequently can be described by a set of conditions. PA1 (1) conservation of mass, PA1 (2) conservation of momentum, PA1 (3) conservation of energy (first law of thermodynamics), and PA1 (4) thermodynamic equilibrium. PA1 (1) Explicit procedures use mobilities and capillary pressures computed as functions of saturations at the beginning of a timestep. The saturations are known from the previous timestep calculations. The mobilities and capillary pressures are assumed to maintain the same values during a timestep that they had at the beginning of the timestep. PA1 (2) Implicit procedures use mobility and capillary pressure calculated as functions of saturation at the end of the timestep. The values are not known until calculations for the timestep have been completed. As a result, they must be determined using an iterative process.
From these principles, a complete mathematical model of the reservoir can be developed which is a combination of (i) boundary conditions, (ii) initial conditions to describe conditions at zero time, and (iii) a system of equations that govern the behavior of the fluids in the reservoir. The governing equations are obtained by combining the following physical principles:
The resulting equations are nonlinear partial differential equations that describe the unsteady-state flow of all fluids throughout the reservoir and relate the pressure and saturation changes of the fluids with time throughout the reservoir. Examples of methods for constructing mathematical models of reservoirs are described in the following books: Peaceman, D. W., Fundamentals of Numerical Reservoir Simulation, Elsevier Scientific Publishing Company, Amsterdam, (1977); and Mattax, C. C. and Dalton, R. L., Reservoir Simulation, Monograph Volume 13, Society of Petroleum Engineers (1990).
Solution of the mathematical model requires determining values of the independent parameters that satisfy all the governing equations and boundary conditions. In general, there are two solution choices: analytical or numerical methods. Analytical methods are not used because the equations are highly nonlinear and are practically impossible to solve by today's methods. The most common method for solving the equations uses a discretization process.
The second step of conventional reservoir simulation (step 2 in FIG. 1) involves subdivision of distance and time into definite, specified increments. This discretization process transforms the continuous partial differential equations to a finite dimensional system of algebraic equations (commonly called finite difference equations). Space is discretized by subdividing the reservoir volume into contiguous cells conforming to a predetermined grid pattern, frequently called "gridblocks" or "gridcells." A typical reservoir model can have thousands of cells. These cells will have sides with lengths ranging from a few meters to a few hundred meters.
The life of a reservoir is discretized (or divided) into time increments. A simulator computes changes in each gridcell (flow, pressure, etc.) over each of many timesteps. Typically, conditions are defined only at the beginning and end of a timestep, and nothing is defined at any intermediate time within a timestep. Consequently, conditions within each gridcell may change abruptly from one timestep to the next. Usually, timesteps are chosen to be small enough to limit sizes of these abrupt changes to acceptable limits. The size of the timesteps depends on accuracy considerations and stability constraints; generally the smaller the timestep the more accurate the solution, but smaller timesteps require more computational work.
Although discretization of the reservoir is an abstraction, it is qualitatively correct to visualize gridcells as well-stirred tanks with permeable sides. The contents of a gridcell are considered uniformly distributed within the cell and the rates at which fluids flow in or out are determined by the permeabilities of the sides of the cell and the pressure differences between adjacent gridcells. In essence, the mathematical problem is reduced to a calculation of flow between adjacent gridcells.
Because of the discretization process, discontinuities in saturation and pressure exist at interfaces between gridcells. Similarly, saturation and pressure are also discontinuous in time, with the reservoir condition defined only at the end of each timestep. Properties that are functions of saturation (such as relative permeability and capillary pressure) can change significantly over a timestep. In the real reservoir, at any particular time, a gridcell has only a single value of each phase saturation and any property that is saturation dependent. To represent variations in reservoir properties, the properties of the gridcells must differ from cell to cell. The abruptness with which a property changes between neighboring cells is largely a function of gridcell size. During a timestep, the saturation-dependent properties will also have assumed many values. However, in reservoir simulation a single value for each saturation-dependent variable must be used to describe conditions over the entire timestep. It is therefore apparent that the finite difference equations can only approximate the partial differential equations that hold for any point within the grid system at a predetermined timestep.
The third step in conventional reservoir simulation (step 3 in FIG. 1) is to linearize the nonlinear terms that appear in the nonlinear finite difference equations and, based on this linearization, construct a linear set of algebraic equations. A linearization represents a term as a straight line function of the quantity or quantities upon which it depends.
The fourth step in conventional reservoir simulation (step 4 in FIG. 1) is to solve the algebraic equations. It is necessary to solve the equations for the unknown values, which in a reservoir containing oil and water will normally include two of the following four variables: oil pressure, water pressure, oil saturation, and water saturation. Other quantities can be derived from these variables, such as oil production rate and water production rate.
Many simulation methods have been proposed. The method chosen can affect the stability and accuracy of the solution and some methods require more computational work than other methods on a per-timestep basis. The methods differ primarily on how they treat the way the reservoir variables (such as pressure and saturation) vary in time. Most methods involve variations of the following two procedures:
A brief overview of six basic methods for solving finite difference equations is given below. These are the Fully Explicit method, the Implicit Pressure Explicit Saturation method (IMPES), the Fully Implicit method, the Sequential Implicit method (SEQ), the Adaptive Implicit method (AIM), and the Cascade method.
Fully Explicit Method
The Fully Explicit method is the simplest solution method because for the purposes of computing flow between cells it assumes that pressures and saturations remain constant at their initial values throughout each timestep. In effect, the inter-cell flows for each timestep are computed without taking into account the way pressures and saturation are changing. As a result, this method is prone to instability. For example, the procedure may predict that the total mass of a fluid within a cell becomes negative. As the timestep increases, the probability of instability increases until, at some sufficiently long timestep, unphysical predictions are inevitable. In reservoir simulation, these considerations limit the timestep to such a small value that the Fully Explicit method is of no practical use.
IMPES Method
In the IMPES method, saturations are assumed to remain constant at their initial values throughout the timestep. Pressures, on the other hand, are assumed to reach their final values very quickly and to remain there throughout the timestep. In this method, the flow rates and the pressures at the end of the timestep are interdependent and must be determined by iteration. The term "implicit" is used because of the interdependence of the flow rate during a timestep and the pressure solution at the end of the timestep. The pressure- and saturation-dependent terms lag behind the pressure terms. The basic procedure is to obtain a single pressure equation by a combination of flow equations. After the pressure has been advanced in time, the saturations are updated explicitly. When the saturations are calculated, new capillary pressures are calculated, which are explicitly used at the next timestep. Implied in this approach is the belief (or hope) that the saturation and pressure values are not changing very rapidly. In some models the method produces good results. However, the method becomes unstable when there is rapid fluid movement in the reservoir (around the wells, for example) or when capillary pressures affect flow strongly. For this reason, IMPES is often not practical.
Fully Implicit Method
The Fully Implicit method, which is more complex than the Fully Explicit and IMPES methods, is unconditionally stable because it treats both pressure and saturations implicitly. Flow rates are computed using phase pressures and saturations at the end of each timestep. In this method, saturations cannot fall below zero because a fluid can flow only if it is mobile at the at the end of a timestep, and fluids are mobile only for saturations greater than zero. Saturations cannot exceed unity because this would imply negative saturations in other phases. The calculation of flow rates and pressure and saturation solutions involves the solution of nonlinear equations using a suitable iterative technique. As the pressures and saturations are solved, the updating of these terms continues using new values of pressure and saturation. The iteration process terminates when the convergence criteria are satisfied.
The Fully Implicit method is stable for all timestep lengths, and guarantees that saturations remain within physical bounds. The unconditional stability of the Fully Implicit method makes it possible for the simulation to take much longer timesteps than would be possible for a non-implicit computation. Although the Fully Implicit method requires an order of magnitude more computing time per timestep than the IMPES method, the timestep can be several orders of magnitude longer, resulting in a net decrease in computing time.
The main drawback of the Fully Implicit method is the amount of computer time that it requires. In terms of computing cost, the method is generally satisfactory in models of single wells or parts of a reservoir, but it can be quite expensive to use in models of entire reservoirs. Several attempts have been made to reduce the computations required, possibly at the cost of accepting a method that does not permit the timestep sizes of the Fully Implicit method. The sequential implicit method, the adaptive implicit method, and the Cascade method have been proposed as ways of reducing the computational time.
Sequential Implicit Method
The Sequential Implicit method (SEQ method) incorporates implicit treatment of saturations, but without solving simultaneously for pressures and saturations. The SEQ method consists of two steps. The first step solves a set of pressure equations in exactly the same way as is done in the IMPES method. This set comprises a single equation per cell, and solving it yields a complete, new pressure distribution at the end of a timestep. In a second step, the pressure distribution is used to calculate the sum of the velocities of all phases at each boundary between cells. These total velocities are used in constructing a set of saturation equations. This set comprises two equations per cell in three phase cases and one equation per cell in two phase cases and is solved simultaneously to yield saturations at the new time. The second step is an implicit solution for saturations using linearized implicit velocities. Saturations in each cell are determined by using implicit (end-of-time-step) values of relative permeabilities in inter-cell fluid flow terms. This method requires simultaneous solution of all saturation equations. The per-timestep computational cost of the SEQ method, though 50 to 100 percent greater than that of the IMPES method, is typically about one-fifth that of the Fully Implicit method. The SEQ method can take much longer timesteps than the IMPES method, and in some models its timestep size approaches that of the Fully Implicit method; in other models, the SEQ method requires up to several times as many timesteps as the Fully Implicit method. Therefore, in some models the SEQ method takes less computing time overall than the Fully Implicit method; in other models, it takes more.
Adaptive Implicit Method
The Adaptive Implicit method (AIM) is sometimes used in an attempt to combine the superior timestep size of the Fully Implicit method with the low computational cost of IMPES. It is based on the observation that only small regions in the typical model require the stability of Fully Implicit computations, while a simpler method would be adequate in the rest of the model. AIM performs Fully Implicit computations only in the parts of the model where they are believed to be needed and uses IMPES in the rest of the model. Typically, Fully Implicit computations are made in only a small part of the model, containing perhaps ten percent of the model's cells.
The discussion above assumes that the Fully Implicit method can take timesteps of arbitrary length. In theory, this is true. Practically, there are difficulties arising from the need to solve a nonlinear system of equations. The equations are solved using an iterative technique called "the Newton method," which makes successive estimates of the solution. This procedure is often referred to as a Newton iteration. When the Newton method works well, each new estimate is closer to the correct solution than its predecessors. Sometimes, it does not work well, and convergence to the correct solution is obtained slowly or not at all. When convergence is poor, the convergence problem is typically confined to a very few cells. Despite this, the fill Newton iteration computation is normally carried out at all cells (or all implicit cells if AIM is in use). A possible way to improve computational efficiency is to, in some manner, focus the computational effort on the few cells having difficulty.
Cascade Method
The Cascade method provides a way of improving the computational efficiency. This method computes saturation one phase at a time, one cell at a time. At a given point in the computation, the cell at which saturation is currently being computed is selected to be one for which all incoming flows of the selected phase have already been computed. Thus, it is only necessary to compute flows of the phase out of the cell. Assuming that upstream weighting is being used, these flows depend on saturation in the current cell only.
The first step in the Cascade method is to compute pressure, just like the IMPES method. Once all pressures are known, it is always possible to find a cell for which all incoming flows of a selected phase have already been computed. Initially, this is likely to be a cell into which the fluid phase is being injected. Once saturation in this cell has been computed, all flows into at least one of its neighbors will be known. Saturation in this cell can be computed, after which all flows into at least one more cell will be known, etc. In practice, the computation is performed in effect by sorting the cells such that the one with the largest potential is first, the one with the next largest potential second, and so on down to the one with the smallest potential, which is last. Then the cell-by-cell computations are performed in this order. Since all flow into a given cell must come from cells having larger potentials, all flows into the current cell are known at any stage of the computation.
The Cascade method achieves its effectiveness by localizing the nonlinear computations to a single cell and, typically, a single unknown within that cell. Because this problem is nonlinear, it must be solved iteratively. Because the computation involves only one unknown, an iterative procedure to solve it can be selected for reliability rather than efficiency. Only cells at which saturation is changing significantly will require more than one iteration, and only cells at which it is changing rapidly will require more than two or three iterations. As a result, the computational effort is concentrated in the cells in which it is needed most.
The Cascade method is particularly suited for models in which flow of all phases is in the same direction. This condition almost never occurs in practical reservoir models, since phase density differences induce countercurrent flow. To some extent, this problem can be overcome by performing a separate Cascade method computation for each phase. However, the method has other drawbacks. First, in the presence of capillary forces, the assumption that flow of a phase out of a cell depends only on that phase's saturation in the cell is no longer valid. Second, an implicit assumption of the Cascade method is that saturation of a given phase is affected only by flow of that phase. In the presence of interphase mass transfer, this assumption is not satisfied. Perhaps because of these shortcomings, the Cascade method is not widely used.
Linearizations
The linearization step (step 3 of FIG. 1) and the solution step (step 4 of FIG. 1) are dependent on each other. In the process of linearization, the algebraic equations will have different forms depending on the solution technique chosen. For example, the IMPES method linearizes only the pressure-dependent terms, such as specific volume. The specific volume is therefore expressed as a linear function of pressure. The SEQ method linearizes the same pressure-dependent terms with respect to pressure, and it also linearizes phase fractional flow terms with respect to saturations. The Fully Implicit method linearizes the pressure-dependent terms with respect to pressure and the saturation-dependent terms (which comprise the relative permeabilities and capillary pressures) with respect to saturations.
The Fully Implicit method uses the linearization differently from the IMPES method and the SEQ method. In the IMPES and SEQ methods, solving the linearized equations gives a solution at the end of each timestep. In the Fully Implicit method, which uses the Newton method, solving the linearized equations yields approximate solutions. Newton iterations are repeated until the resulting estimates of the solutions are considered to be accurate enough based on pre-specified convergence criteria.
One difficulty with conventional linearization procedures is that under certain reservoir conditions the computed saturations can be negative, a physical impossibility. If subsequent calculations are made based on the impossible saturations, erroneous results are produced. These erroneous results can occur, for example, in a gridcell having a high fluid throughput during a timestep. This occurs frequently near wells, or in gridcells that simulate fractures in a reservoir, where gridcells are small and the fluid velocities are large. Hundreds or thousands of gridcell pore volumes can flow through a gridcell in a single timestep. This high throughput can cause stability problems in reservoir simulation.
The most common way to deal with this potential instability is to use the Fully Implicit method, using a Newton iteration at each timestep with damping used to prevent the iteration from diverging. Because a coupled set of equations is being solved, the Fully Implicit method can be expensive computationally. The SEQ method has also been proposed for high throughput gridcells, but it is typically less stable and usually takes smaller timesteps than the Fully Implicit method.
A need exists for an improved reservoir simulation computation method that is more reliable and accurate than the computational methods used in the past and at the same time is computationally efficient.